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We invoke a metric to quantify the correlation between any two earthquakes. This provides a 
simple and straightforward alternative to using space-time windows to detect aftershock sequences 
and obviates the need to distinguish main shocks from aftershocks. Directed networks of earthquakes 
are constructed by placing a link, directed from the past to the future, between pairs of events that 
are strongly correlated. Each link has a weight giving the relative strength of correlation such 
that the sum over the incoming links to any node equals unity for aftershocks, or zero if the event 
had no correlated predecessors. A correlation threshold is set to drastically reduce the size of the 
data set without losing significant information. Events can be aftershocks of many previous events, 
and also generate many aftershocks. The probability distribution for the number of incoming and 
outgoing links are both scale free, and the networks are highly clustered. The Omori law holds 
for aftershock rates up to a decorrelation time that scales with the magnitude, m, of the initiating 
shock as tcutoft ~ 10''™ with /3 ~ 3/4. Another scaling law relates distances between earthquakes 
and their aftershocks to the magnitude of the initiating shock. Our results are inconsistent with the 
hypothesis of finite aftershock zones. We also find evidence that seismicity is dominantly triggered 
by small earthquakes. Our approach, using concepts from the modern theory of complex networks, 
together with a metric to estimate correlations, opens up new avenues of research, as well as new 
tools to understand seismicity. 



I. INTRODUCTION 



Seismicity is an exceedingly intermittent phenomena 
|22'| exhibiting strong correlations in space and time. 
Since the seismic rate increases sharply after a large 
earthquake in the region, events typically are classified 
as aftershocks or main shocks, and the statistics of af- 
tershock sequences are studied. Usually, aftershocks are 
collected by counting all events within a fixed space-time 
window 0, 1271 following a main event. The 
size of the window may vary or scale with the magnitude 
of the main event |23l |. and other refinements have been 
made ^sj- However, this method does not allow a like- 
lyhood to be estimated that an event thereby collected 
is actually correlated to the main event under consid- 
eration. As a result, no straightforward algorithm ex- 
ists to decide if the space-time windows are too large or 
too small for minimizing errors in the procedure. Also, 
the method cannot be easily extended to examine remote 
triggering pH Il7| where main shocks may trigger after- 
shocks at great distance in space, or perhaps far away in 
time. In addition, aftershocks may have several preceding 
events to which they are correlated, perhaps with over- 
lapping space-time windows. Using conventional meth- 
ods the conjecture that aftershocks can rumble on for 
centuries p3j cannot be tested. 

In fact, a growing body of work indicates that the 
distinction between aftershocks and main shocks is rela- 
tive 0, 01 ■ There is no unique operational way to distin- 
guish between aftershocks and main shocks p. They are 



not caused by different relaxation mechanisms IT^ . 
Besides, a strict distinction may not be the most useful 
way to describe the dynamics of seismicity. 

A particular nuisance with employing space-time win- 
dows arises from the entanglement of a vast range of 
scales in space, time and magnitude in seismici ty. One 
scale-free property is the Gutenberg and Richter [l^l (G- 
R) distribution for the number of earthquakes of magni- 
tude m in a seismic region, 
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with 6 usually « 1. A second is the fractal appearance 
of earthquake epicenters 0, |22, 113 , where the fractal 
dimension dt « 1.6 in S. California % A third is the 
Omori law [33,|39| for the rate of aftershocks in time. 
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where c and K are constant in time, but depend on the 
magnitude m |3^ of the earthquake. 

One way forward has been suggested by who take 
the perspective of statistical physics: Neglecting any clas- 
sification of earthquakes as main shocks, foreshocks or af- 
tershocks, analyze seismicity patterns irrespective of tec- 
tonic features and place all events on the same footing. 
They consider spatial areas and their subdivision into 
square cells of length L. For each of these cells, only 
events above a threshold magnitude m are included in 
the analysis. In this way, one can obtain a distribution 
of waiting times 0,0,111 and distances between succes- 
sive events with epicenters both in the same cell of linear 
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extent L. Since both the threshold magnitude and the 
length scale of the cell (or the space window) are arbi- 
trary, one looks for robust or universal features of this 
distribution that appear when these parameters are var- 
ied. 

We have previously introduced an alternative 
method for characterizing seismicity that completely 
avoids using fixed space-time windows and, at the same 
time, makes available powerful concepts and methods 
that are being developed to describe complex networks 
QilEH^' Here we extend the method to allow more than 
one incoming link per earthquake. This allows the net- 
work to deviate from a tree structure and exhibit charac- 
teristic properties of complex networks, such as "cluster- 
ing" 3] , which may be relevant to characterizing seismic- 
ity. Our method also takes into account, in an unbiased 
way, that an aftershock can be correlated to many pre- 
vious events. By unbiased, we mean that we do not fix 
any length, time or magnitude scales for identifying af- 
tershocks. Nor do we fix the number of events they can 
be aftershocks of. 



II. THE METHOD 

A general description of our method is as follows: The 
first step is to propose, as a null hypothesis 20], that 
earthquakes are uncorrelated in time. Then we detect 
instances when that hypothesis is strongly violated, in- 
dicating that the opposite is true. The second step is 
to assign a real number, or metric, that quantifies the 
correlation between any two earthquakes, based on gross 
violations of the null hypothesis. The third step is to 
construct a directed network where the events that are 
correlated according to the metric are nodes connected 
by links. Each link contains several variables such as the 
time between the linked events, the spatial distance be- 
tween their epicenter or hypocenters, the magnitudes of 
the earthquakes, and the metric or correlation between 
the linked pairs. We can study the statistical properties 
of the network and its ensemble of space/time/magnitude 
variables to gain new insights into seismicity. Note that 
many variations of the null hypothesis and associated 
metric are possible, but the key feature of a useful null 
hypothesis, in this context, is that earthquakes are un- 
correlated in time. 

The null hypothesis, which we previously used, is 
that earthquakes occur with a distribution of magnitudes 
given by the G-R law, with epicenters located on a frac- 
tal of dimension dj, at random in time. Of course, it is 
patently false that earthquakes are uncorrelated in time. 
It is also unclear if epicenters form a monofractal with di- 
mension df < 2. The point is to look for strong violations 
of the null hypothesis. 

Consider an earthquake j in the seismic region, which 
occurs at time Tj at location Rj . Look backward in time 
to the appearance of earthquake i of magnitude rui at 
time Ti, at location i?;. One can ask, how likely is event 



i given that event j occurred where and when it did? 
According to the null hypothesis, the expected number 
of earthquakes of magnitude within an interval Am of 
rui that would be expected to have occurred within the 
time interval t — Tj — Ti seconds, and within a distance 
I — I i?i — Rj I meters is 

n,j = [const) tl"^! io~^™* Am . (3) 

Note that the space-time domain (i, I) appearing in Eq.|3| 
is selected by the particular history of seismic activity in 
the region and not preordained by any observer. The 
constant term in Eq. |31 is estimated by the overall seis- 
mic rate in the region over the time span of recorded 
events and is evaluated later. However, our results are 
insensitive to the precise value of this constant, since its 
value is absorbed into a threshold we define later, c<. 
We find that many of the statistical properties of the 
networks are robust with respect to varying parameters 
such as c<, df, and h. In particular, we can choose dt_= 2 
without substantially varying the results. See also 

Consider a pair of earthquakes where <C 1; so 
that the expected number of earthquakes according to the 
null hypothesis is very small. However, event i actually 
occurred relative to j, which, according to the metric, is 
surprising. Hence, it is unlikely that the pair would oc- 
cur in that space-time domain if they were uncorrelated. 
A small value ^ 1 indicates that the correlation be- 
tween j and i is very strong, and vice versa. By this ar- 
gument, the correlation Cij between any two earthquakes 
i and j can be estimated to be inversely proportional to 
n,y , or 

C^, = I/Uij . (4) 

As we show later, the distribution of the correlation vari- 
ables Cij for all pairs i,j is extremely broad. Therefore, 
for each earthquake j , some exceptional events in its past 
have much stronger correlation than all the others com- 
bined. These strongly correlated pairs of events can be 
marked as linked nodes, and the collection of linked nodes 
forms a sparse network of highly clustered graphs. Unless 
otherwise stated in this work, earthquakes are linked only 
if their correlation value, Cy, is greater than c< = 10*, 
or the expected number of events according to the null 
hypothesis, riij, is less than lO"**. The error made in 
ignoring weakly linked pairs of events is discussed later. 

In the language of modern complex network theory 
QIEHOli a time-oriented weighted network grows, where 
nodes (earthquakes) have internal variables (magnitude, 
occurrence time, and location), and links between the 
nodes carry a strength (the correlation Cij) and are di- 
rected from the older to the newer nodes. Empirically, we 
find that both the distribution of outgoing and incoming 
links are scale free. The network is composed of highly 
clustered, disconnected graphs of correlated earthquakes. 
Events with incoming links, or aftershocks, typically con- 
nect to many previous events rather than just one. How- 
ever, the networks are sparse and the number of links in 
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the network is much less (about 0.1%) than the number 
of pairs of earthquakes. We find neither that every earth- 
quake is correlated to every other event, nor that events 
typically are correlated to zero or one previous events, 
but a picture in between where the number of events an 
aftershock is correlated to is scale-free. 

Due to the continuous nature of the link variable, c^-, 
no event is purely an aftershock or a main shock, and it is 
not possible to separate events into distinct classes. This 
is consistent with previous studies indicating no physical 
distinction between main shocks and aftershocks [a, ■ 
Note that singularities in Eq. 13 are eliminated by taking 
a small scale cutoff in time (here tmin = 60 sec) and a 
minimum spatial resolution (here Zmin = 100 meters). 



III. DATA AND PARAMETERS 

The catalog we have analyzed is maintained 
by the Southern California Earthquake Data Cen- 
ter, and can be downloaded via the Internet at 
[http://www.data.scec.org/ftp/catalogs/SCSN/ We use 
data ranging from January 1, 1984 to December 31, 2003, 
and follow a procedure similar to our previous work (see 
Baiesi and Paczuski 0| for more details). 

The relevant quantities for our present work are sum- 
marized in Tables ^ El and IIIII Events with magnitude 
smaller than m< — 3 are discarded, and Am — 0.1. 
The number of earthquakes or nodes in the network con- 
structed using the entire catalog is A'nodo = 8858. The 
&- value of the G-R law is 6 ~ 0.95 for this data set, while 
d/ ~ 1.6 was found by Corral f3|. 

We consider two closely related variants of the met- 
ric: in the two-dimensional (2D) version, the earthquake 
depth is not considered, and the distance between two 
events i and j is measured as the arc length on the Earth's 
surface, 
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where the Earth radius is Rq — 6.3673 x 10^ meters, and 
{6i,(j)i) are the latitude and longitude, in radians, of the 
epicenter of th z'th event in the catalogue. 



TABLE L Network quantities for the Southern California data 
set, unless otherwise noted. 



The second version (3D) takes into account the depth 
hi of each event. Hence Euclidean distances between 
hypocenters are calculated. 



hj — 



\ 



„a\2 



(6) 



with 
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and df in the metric is replaced by the hypocenter fractal 
dimension Df, which is approximated asl?/=ci/-|-l~ 
2.6 for Southern California. Thus, the 3D metric is 
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Most of the statistical results we find are not sensitive to 
the choice of the metric, nor to the precise values of b, 
df, or Df. For this reason, we pick as a standard metric 
the 2D version (Eq. ^ , and use the 3D metric only when 
explicitly stated. 

The constant in Eq. |31 was estimated to be const = 
10^^^ for the 2D metric using the same method as in 
Baiesi and Paczuski P|. However in that work a fractal 
dimension of d/ = 1.2 was inadvertently used to estimate 
const, resulting in a different value. Similarly, here we 
compute const' — 10~^^ for the 3D metric, Eq. ©. Both 
values give consistent results, but they are not expected 
to be precise due to the high variability of seismicity rates 



TABLE IL Parameters for the network of Southern California 
obtained with the 2D version of the metric, unless otherwise 
noted. 



Quantity 



Symbol Value 



seismicity constant (see Eq.|3 const 10~^^ 

correlation threshold c< lO* 

number of links Mink 166507 

average in-degree (fcin) 18.8 

number of clusters A'^ciustcr 2252 

TABLE IIL Parameters for the network of Southern Califor- 
nia obtained with the 3D version of the metric, unless other- 
wise noted. 



Quantity Svmbol Value 


magnitude threshold m< 


3 


Quantity 


Symbol Value 


magnitude precision Am 


0.1 


seismicity constant (see Eq.|SJ 


const' 


10-15 


Gutenberg-Richter exponent b 


0.95 


correlation threshold 


c< 


10* 


fractal dimension of epicenters df 


1.6 


number of links 


Mink 


154792 


fractal dimension of hypocenters Df 


2.6 


average in-degree 


(fcin> 


17.5 


number of earthquakes A'^node 


8858 


number of clusters 


cluster 


2327 
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in the region even over a time span of years. However, 
varying the constants, const and const ', in our analysis is 
equivalent to varying the correlation threshold for link- 
ing events, c<. We observe that many of the statistical 
results presented below are robust to variations of c<. 
This is primarily for two reasons. First, as we will show 
the distribution of link weights P{c) is very broad and 
doesn't pick out preferred values. Second, the distribu- 
tions we compute are weighted using the link weight. Re- 
ducing the threshold for included links only adds earth- 
quake pairs that give progressively lower contribution to 
the final correlation structure. We give later a numer- 
ical estimate for the error made in throwing out these 
degrees of freedom. The advantage, obviously, is that a 
sparse network with c< chosen appropriatly, enables us 
to vastly reduce the size of the data set from approxi- 
mately 10-100 Gigabytes to around 10 Megabytes or so, 
without losing important information about correlations 
between earthquakes. The network allows a 'rcnormal- 
ization' which removes irrelevant degrees of freedom, or 
links with low weights while keeping important ones. 



IV. RESULTS 

Networks constructed using our method are shown in 
Fig. n For comparison, networks obtained with the 
2D metric [Fig. QJa) and (b)] and with the 3D metric 
[Fig. ^c) and (d)] are both displayed. For visual clarity, 
a higher threshold for earthquake magnitudes to< = 4 
was used in order to reduce the number of nodes and 
links in the figure. Adjusting the parameter c< slightly, 
two networks with a similar number of links, with very 
similar clusters of correlated events are formed. There is 
a more abundant presence of long distance links in the 2D 
version but the similar details in the Northridge clusters 
[Fig. njb) and (d)] suggest that it is mainly seismic his- 
tory that determines the network structure, rather than 
the precise details of our metric. 



A. Explanation of Method 

Fig.|21shows the probability distribution of correlation 
values, -P(c), obtained by sampling the values Cy over 
all earthquake pairs in the data set. It is a fantastically 
broad distribution that exhibits power law behavior over 
sixteen orders of magnitude (in the 3D case) : 

P(c) ~ c-^ (9) 

with T = 1.43 ± 0.03 using the 2D metric and r = 1.38 ± 
0.03 using the 3D one. 

Given such a broad distribution, for any earthquake 
J, some extreme events i exist whose correlation Cij are 
much larger than all the others. Therefore, it makes sense 
to represent these earthquake pairs as nodes that are 
linked, while not linking pairs that have much smaller 



values of Cij. Then the sequence of earthquakes may be 
usefully represented as a sparse network, where links ex- 
ist between the most strongly correlated events, i.e. those 
pairs (i,j) where Cij > c<. Hence a natural decomposi- 
tion of the network into disconnected clusters is achieved, 
where the first earthquake in the directed cluster has no 
incoming link, or correlation variable into it greater than 
c<. Clearly, the first earthquake in the entire catalogue 
also no incoming link. The correlated events are reliably 
detected when c< is greater than one but not extremely 
large. In the latter case, correlated events detach, and a 
very fragmented network appears. For small c< some un- 
correlated events make links, and a giant cluster appears. 
Both for the 2D and for the 3D case we set c< = 10*, un- 
less otherwise noted, obtaining a similar number of links 
in the realization of the networks. 



B. The Scale- Free Network 

The resulting network of earthquakes is scale free. As 
shown in Fig. |51 both the distribution of the number of 
incoming links, or the "in-degree" fcin, to a node and the 
distribution of the number of outgoing links, or the "out- 
degree" kout, to any node exhibit power law behavior, 

P(fci„) ~ l/kn , P(fcout) 1/fcout (10) 

up to a degree « 100. 

1. Aftershocks with more than one main shock 

Since an earthquake can have more than one incoming 
link, in attributing aftershocks to an event we must be 
careful not to overweight aftershocks with many incom- 
ing links. To prevent the overcounting of aftershocks, 
one can consider a new event with two incoming links, 
for example, to be "half an aftershock" of both of its 
precursors, or they can be weighted in a different fashion 
according to their correlation values. In general, we can 
attribute the relative correlation to previous events, so 
that each event contributes a total weight of unity to the 
global aftershock number if it is linked to at least one 
previous event, and zero otherwise, as follows: 

For each event j that has at least one incoming link, 
so that it can be called an aftershock, define a weight for 
each "parent" earthquake i it is linked to as 



where the sum is over all earthquakes k with links going 
into j. The weighted number of aftershocks attributable 
any event i is then 

out 

iVaftor,i = ^ Wij . (12) 
J 
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FIG. 1: (a) Scale free network of earthquakes obtained with the 2D metric using m< — 4 and c< = 10''. This network has 
791 nodes, and its 7931 links follow the color code in the legend. Several clusters are evident, the biggest being related to the 
Landers earthquake. The Northridge cluster, enclosed within a dashed box in (a), is enlarged in (b), where the solid black 
dot represents the epicenter of the Northridge event. Figures (c) and (d) are obtained using the 3D metric with m< = 4 and 
c< = 5 X 10^ giving 7947 links, very close to the number of links in the 2D version. Note that the networks found using these 
two metrics are similar, indicating that the method is robust to variations in the metric. 



Here, the sum is over all of the outgoing links from event 
i. In the limit that rj —>■ oo, the extremal network stud- 
ied by Baiesi and Paczuski [J is recovered, since only the 
single incoming link to aftershock j, with the largest cor- 
relation Cij, contributes. In that case, for each node, the 
quantity A'^aftcr discussed here coincides with the quan- 
tity fcout in Baiesi and Paczuski In the following, we 
consider the case rj = 1. 

Sampling over all earthquakes, we get a probability 
distribution for the number of (weighted) aftershocks 
as also shown in Fig. |21 It is a power law distribu- 
tion, P{Naftci) (^aftcr)"''' Scaling over more than three 
decades, with an index 7 = 2.0(1). This distribution is 



very close to the distribution we obtained previously for 
the number of aftershocks in the extremal network, cor- 
responding to the limit rj ~> 00. The distribution for the 
number of weighted aftershocks appears to be universal, 
in the sense that the power law exponent does not depend 
on 77, for 77 > 1. Note that chosing a positive rj gives more 
weight to more strongly correlated pairs and is therefore 
consistent with using a threshold c< to eliminate weakly 
correlated ones. 
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FIG. 2: The probability distribution of the correlation, c, 
between all earthquake pairs in the data base, with m< = 3, 
using both the 2D metric and the 3D one. They are scale free 
distributions over many orders of magnitude. The threshold 
c< = 10* where correlations are considered significant and 
links are made is indicated in the figure. Note that, with that 
threshold, most links are eliminated from the network, giving 
a reduced data set to examine seismic properties. 



Or 



a. 



o 


h 

in 


♦ 


k , 

out 


▲ 





♦ 

o 



-6- 




log,n degree 



'■^ 

3 



FIG. 3: The in-degree and out-degree distributions of the net- 
work of earthquakes and aftershocks. The out-degree, fcout, is 
the number of outgoing links from an earthquake, linking it to 
its aftershocks. The in-degree, fcin, is the number of incoming 
links to an earthquake, linking it to its main shocks. These 
two distributions are similar. Also shown is the distribution 
of the weighted number of aftershocks, P(A'aftcr), from any 
event, using 77 = 1 in Eqs. 1111 and 1121 This has the same 
scaling behavior as the extremal network, with 77 00. 



2. Is Seismicity driven by small or large earthquakes? 

The average number of aftershocks of an earthquake of 
magnitude m has been proposed to scale with m |^ [s^ 
as 



^aftor(m) - 10° 



(13) 



Larger main shocks release more energy and therefore 
"trigger" more aftershocks than smaller earthquakes. 
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FIG. 4: The total number of outgoing links from all events of 
magnitude m (diamonds), the total number of weighted after- 
shocks of all events of magnitude m (triangles) and average 
number of incoming links ((fcin(m)), circles and histogram) 
as a function of m. The first two quantities have a scaling 
consistent with the law ~ 10"™, with a' ^ b - 0.5 ^ 0.45 
for outgoing links, and a ~ b ~ 0.15 ~ 0.8 for weighted after- 
shocks. For (fein(m)) there is no evident departure from the 
average value (hn) ~ 18.8 (dotted line). 



However, smaller earthquakes are more frequent, as in- 
dicated by the G-R relation. The total number of after- 
shocks generated by all earthquakes of magnitude m is 
therefore given by the product of these two relations as 



WTOT(m) = TVaftcr (m)P(m) ^ lO^"-")™ 



(14) 



If the exponent a > b then small earthquakes are the 
dominant triggering mechanism for seismicity, whereas if 
a < b the large earthquakes dominate aftershock produc- 
tion. Often it is assumed that a — b ^21. J6J . 

Recently, Helmstetter analyzed earthquake cata- 
logues by means of a "stacking" method using space-time 
windows, and found that aftershocks were predominantly 
triggered by small earthquakes. She determined the value 
of the exponent a to be between a = 0.72 and a = 0.82, 
depending on the parameters of the aftershock detection 
algorithm that she used. 

Fig. 0] shows A^TOT(m-) obtained using the 2D met- 
ric with r; = 1. The results shown in this figure also 
suggests that small earthquakes are the dominant mech- 
anism driving aftershock production. We determine the 
value of the exponent a w 0.8, consistent with Helmstet- 
ter's previous findings. Thus at this rather detailed level, 
results obtained using our method are consistent with 
results obtained using traditional methods of aftershock 
detection. This is true despite the fact that aftershocks 
in our algorithm are typically attached to many previous 
events rather than just one, and no space-time scales are 
used by us for aftershock identification. 

Fig. ^ also shows the total number of links emanating 
from events of magnitude to, A^out('7i). This corresponds 
to an unweighted aftershock number, with A^out('7i) ~ 
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lQ-{a -b) gj^^ Q,/ _ g_45^ Since both a and a' are less 
than 6, our results suggest that irrespective of the man- 
ner in which aftershocks are weighted, small earthquakes 
are the dominant mechanism driving aftershock produc- 
tion. Note however, that the largest events may appear 
to present a deviation from this behavior, but the statis- 
tical uncertainties of single events are large. 



3. Bath's Law 

In Fig. 21 we also show the dependence of the number 
of incoming links to a node on its magnitude. The quan- 
tity (fcin(m)), is the average number of incoming links 
to earthquakes of magnitude m. This quantity is inde- 
pendent of earthquake magnitude for 3 < m < 5. For 
larger magnitudes, the poor statistics forces us to av- 
erage fci,i(m) over wider bins, chosen so that there is a 
significant number of events inside each one. These av- 
erages are indicated as horizontal lines in Fig. 21 while 
vertical lines denote the bin boundaries. The averages 
do not show any detectable trend for larger earthquakes. 

Our results suggest that earthquakes of all magnitudes 
are equally likely to be aftershocks, and support the con- 
clusion reached by Helmstetter and Sornette that 
observations of Bath's law are due to biases in labellini 
earthquakes as aftershocks. According to Bath's law 
the average magnitude difference between a main shock 
and its largest aftershock is around 1.2, independently 
of the main shock magnitude. Of course, the definitions 
we use here for main shocks, as nodes giving outgoing 
links, and aftershocks, where nodes have incoming links 
(so that a single event can be both a main shock and 
aftershock), differs from the standard definition. 



4- Clustering of nodes 

Among the concepts in network theory that may be 
useful to characterize seismicity, and are not accessible 
via other approaches, the clustering of nodes deserves 
particular attention. Indeed, the clustering in space and 
time of earthquakes can be quantified in their network 
by the clustering coefficient. The clustering coefficient of 
a node i is the number of linked triangles it forms 
with its ki neighbors (or equivalently, the pairs of linked 
neighbors) divided by the maximum number of linked 
triangles it could potentially have {ki{ki — l)/2), i.e., 



2A,; 



ki(ki 1) 



(15) 



This definition ignores the directionality of links. Thus, 
in this formula the degree of node i, is the sum of its 
incoming and outgoing degrees, i.e. ki = fci,in + ^i.out- In 
all cases < Ci < 1, and = if less than two links 
are joined to node i, or if no links between its linked 
neighbors are present, while Ci — I only if all neighbors 
are linked to each other. 
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FIG. 5: The network density, or clustering coefficient as a 
function of the degree, k, of a node, where k — kin + kout- 
For small values of k, the clustering coefficient is independent 
of k and c<, and C ~ 0.80. For large values of k, C{k) tends 
toward power law behavior C{k) ~ k'^ with 5 « 1.3. The 
power law regime takes place at smaller k for larger thresholds 
c<. 



Using Eq. (|15|l to compute the average clustering co- 
efficient of the network. 



c 



(16) 



i=i 



we obtain C = 0.50 for to< = 3. This value is relatively 
stable with respect to variations of r7i<. For instance, 
with m< = 4.5 we get C = 0.55. The same values are 
obtained for the 3D version of the metric. These are 
remarkably high values of C, compared to many other 
complex networks, such as technological or biological 
ones [30j |. 

a. Universal clustering properties of seismic networks 
The average of the clustering coefficient can be performed 
over nodes with the same degree k. This quantity is 
shown in Fig. jsl where one can observe that it does not 
depend on k for small values of k and approaches a power 
law C{k) ^ k^^ with i5 « 1.3 for large values of k. This 
power law behavior is typically found in networks with 
a modular structure |135|- At small k, C{k) approaches a 
universal value approximately equal to 0.8, which is in- 
dependent of k and of the thresholds c< and r7i< used 
to construct the network. The power law exponent 6 
at large values of k also appears to be independent of 
c< , and may also be a universal quantity for seismic net- 
works. 



C. Scaling Law for Aftershock Distances 

In the network constructed using the 2D metric, the 
link length, I, is the distance between the epicenters of 
an earthquake and one of its aftershocks, weighted ac- 
cording to the link weight w. In the corresponding net- 
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work constructed using the 3D metric, the Hnk length is 
the distance between the hypocenters of an earthquake 
and one of its aftershocks, weighted according to the hnk 
weight w. The distribution PmQ-) of hnk lengths depends 
on the magnitude m of the predecessor, being on aver- 
age greater for larger m. To compute this distribution, 
we put the weight of each link into a bin corresponding 
to its I value and the magnitude of the predecessor m 
to get Pm{l)- A maximum in the distribution occurs, 
which shifts to larger I on increasing m, as shown in 
Fig. This behavior is superficially consistent with us- 
ing larger space-time windows to collect aftershocks from 
larger events, or the Kagan (23| hypothesis of aftershock 
zone scaling with main shock magnitude. 



(a) 2D metric 
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1. Comparison with Aftershock Zone Scaling 

It is widely believed that an aftershock zone exists 
which is equivalent to the rupture length. Within the 
aftershock zone, earthquakes generate aftershocks, while 
outside the zone they do not. The rupture length, R, is 
believed to scale as i? ~ iQO.Sm ^[^]^ ^j^g magnitude of the 
main shock. This is a restatement of the relation derived 
by Kanamori and Anderson j^^, who argued that the 
seismic moment M ~ iQ^-5m g^g^les with R as M ^ R^, 
at least for intermediate magnitude earthquakes. For a 
generalization to all earthquakes, see Kagan |2^. In this 
scenario main shocks of all magnitudes generate after- 
shocks at the same rate within their respective aftershock 
zones, so that the greater number of aftershocks coming 
from large events is due solely to their larger aftershock 
zones. Needless to say, the observation of aftershock zone 
scaling is based on the idea that the aftershock zone is 
finite - on the order of tens of kilometers for large main 
shocks. 

In contrast, we find the distribution of lengths between 
main shocks and their aftershocks exhibits no cutoff at 
large distances, but rather decays slowly as a power law 
with Z, up to the linear extent of the seismic region cov- 
ered by the catalog, hundreds of kilometers. The two 
distributions are both consistent with a scaling ansatz: 



(17) 



where I is measured in meters and F(x) is a scaling func- 
tion. Remarkably in both cases, a « 0.37. Note in par- 
ticular that cr 7^ 0.5. 

For x ^ 1, the tail of the scaling function is a power 
law, i.e. F{x) ~ x"^ with A w 2 (2D) or A « 2.6 
(3D). The results obtained using the data collapse tech- 
nique applied using ansatz H17|l are shown in the insets 
of Fig. Ela) and (b). Such slow decays at large dis- 
tances calls into question the use of sharply defined space 
windows for collecting aftershocks, as already pointed 
out by 113. The length scale we find l*^ ~ ioO-37m 
describe the fat-tailed distribution of distances between 
earthquakes and their aftershocks should not be confused 
with the scaling of a finite aftershock zone as proposed 
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(b) 3D metric 
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FIG. 6: Link length distribution for different magnitudes of 
the emitting earthquake, (a) for the 2D case and (b) for the 
3D one. The length where the maximum in the distribution 
occurs increases with magnitude roughly as imax ~ iQO.iVm 
both cases. Both distributions also have a fat tail, extending 
up to hundreds of kilometers even for intermediate magnitude 
events. These distributions are consistent with a hierarchical 
organization of events, where big earthquakes preferentially 
link at long distance with intermediate ones, which in turn 
link to more localized aftershocks, and so on. Insets: Dis- 
tributions rescaled according to Eq. 1171 with a = 0.37 and 
with m equal to the central magnitude of the range for each 
distribution. 



by Kagan 23]. Instead, our results are consistent with 
observations of remote triggering of aftershocks by Hill 
et al. d7j and Gomberg et al. Jjj as well as the obser- 
vation that the distribution of distances between subse- 
quent earthquakes in regions of size L is a power law, not 
trivially given by the correlation dimension, ^2 of earth- 
quakes, and which is cutoff only by the size of the region 

L a. 
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FIG. 7: The Omori law for aftershock rates. These rates are 
measured for aftershocks Unked to earthquakes of different 
magnitudes. For each magnitude, the rate is consistent with 
the original Omori law, Eq.|51 up to a cutoff time that depends 
on m. As guides to the eye, dashed lines represent a decay 
~ The dense curves represent the fits obtained by means 
of Eg. ll9l for m — 3, m = 4, and m — 5. 



D. The Omori Law for Earthquakes of All 
Magnitudes 

Figure [3 shows the rate of aftershocks for the Landers, 
Hector Mine, and Northridge events, obtained with the 
2D metric. The weights, w, of the hnks to aftershocks 
occurring at time t after one of these events are binned 
into geometrically increasing time intervals. The number 
of weighted aftershocks in each bin is then divided by the 
temporal width of the bin to obtain a rate of weighted 
aftershocks per second. The same procedure is applied 
to each remaining event, not aftershocks of these three. 
An average is made for the rate of aftershocks linked to 
events having a magnitude within an interval Am of m. 
Figure[7|also shows the averaged results for m = 3 (1871 
events), m = 4 (175 events), m — 5 (28 events) and 
m = 5.9 (4 events). 

The collection of aftershocks linked to earthquakes 
of all magnitudes is one of the main results of our 
method. Even intermediate magnitude events can have 
aftershocks that persist up to years. Earthquakes of all 
magnitudes have aftershocks which decay according to 
the Omori law [H Hi], 

K 

'^(*) rz ' fo^' ^ < ^cutoff (18) 

c + t 

where c and K are constant in time, but depend on the 
magnitude m 39] of the earthquake. We find that the 
Omori law persists up to a decorrelation time ^cutoff that 
also depends on m. 

A rough estimate of the decorrelation times can be 
extracted by non-linear fits of logj^o ^in{t) vs logj^Q using 
an interpolating function 

Vra{t) ~ t-le-*/*--« . (19) 



FIG. 8: Decorrelation time, fcutoff of the aftershock rates to 
fall out of the Omori regime, as a function of the earthquake 
magnitude. The horizontal line indicates 20 years, i.e. the 
time span of the catalogue. The dashed line is the interpola- 
tion given in Eg. 1201 

The range of the fit excludes short times, where the the 
aftershock rates are not yet scaling as \/t. The short time 
deviation from power law behavior is presumably due to 
saturation of the detection system, which is unable to 
reliably detect events happening at a fast rate. However, 
this problem does not occur at later times, where the 
rates are lower. Some examples of these fits are also 
shown in Fig. [7| for the intermediate magnitude events. 

In Fig.|Slwe show the resulting values of tcutoff, for m 
up to 6. The horizontal dotted line represents the time 
span of the catalogue we study, which precludes accurate 
estimates of much longer tcutoff- Thus, from the rates of 
aftershocks of events with 3 < m < 4.6, where the time 
span of the catalogue is comparable or longer than the 
estimated decorrelation time, we find that the increase of 
^cutoff with TO can be fitted by the function 

teutoff(m) ~ 10^-25+0-74m ^ (2O) 

represented as a dashed line in Fig. [S] It roughly cor- 
responds to tcutoff ~ 11 months for m = 3, and to 
tcutoff ~ 5 years for m = 4. An extrapolation yields 
tcutoff ~ 1400 years for an event with to = 7.3 such as 
the Landers event! However, we stress that Eq.[5n|is just 
rough estimate of tcutoff (™)- 

Note that Helmstetter [iJl also found an Omori law 
for aftershocks of earthquakes of all magnitudes using 
finite space-time windows. However, for this reason, she 
was not able to estimate the decorrelation of aftershocks 
or the cutoff in the duration of the Omori regime for 
different magnitudes. 

V. DISCUSSION 

At present, we are unaware of any reliable method to 
determine the best metric. Thus, the best route to study 
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how sensitive the results are to variations of the metric or 
to the parameters of the metric. Although we have not 
yet made an exhaustive and detailed study to determine 
which properties may be universal and hold for many 
different metrics, several general conclusions are already 
apparent. 



A. Robustness with respect to changes in 
parameters 

Many of the statistical results we find are relatively 
robust with respect to variations in the metric. For in- 
stance, using both the 2D metric (Eq. I^J and the 3D 
metric (Eq. |SJ), similar networks are found as indicated 
qualitatively in Fig. 1. In addition, the scaling behaviors 
demonstrated in Figs. 2-7 are independent of the metric, 
with the notable exception of the exponent A character- 
izing the fat tailed distribution of aftershock distances. 
However, the exponent o for the rescaled variable com- 
bining main shock magnitude and aftershock distance is 
independent of the metric, as is the Omori behavior. Fur- 
thermore, the distribution of correlations Pic) depends 
only weakly on the metric, and the scale-free and clus- 
tering properties of the network are insensitive as well. 

One could object that the values of 6, d/ and D/ can 
depend on the region of the Earth being considered, or 
may fluctuate depending on the specific fault zone being 
studied. However, the statistical results we find, as shown 
in the Figures, are also robust to variations in either of 
these parameters, or of the threshold m<. This robust- 
ness was also found in the Baiesi and Paczuski stud- 
ies of the extremal earthquake network. For instance, 
varying df over a wide range, from 1 to 2 does not al- 
ter considerably the distribution of incoming or outgoing 
links. The distribution of correlations -P(c), is even more 
insensitive to variations of h and d/. Also the Omori law 
with p « 1, shown in Fig. [7| does not depend sensibly 
on the parameters, and holds for aftershocks linked to 
earthquakes of all magnitudes. 

Our interpretation of this observed robustness is that 
the correlation structure of seismicity is unambiguous 
and clear-cut, and has a network structure similar to 
other complex networks. Even if we use an approximate 
measure, or metric, the underlying correlations are suf- 
ficiently strong that they survive the approximation and 
can be reliably detected. 



B. Errors and Data Set Reduction 

One could also object that the parameter c< is ar- 
bitrary, and its choice plays a similar role to choosing 
space-time windows in the traditional manner. How- 
ever, one can consider all pairs of earthquakes, using their 
weights Cij, so the parameter c< is conceptually unnec- 
essary. This differs from the necessity of choosing space- 
time windows in the traditional approach. However, as a 



practical matter, to reduce the size of the data set, it is 
useful to choose a particular c< , and thereby construct a 
sparse network. The choice involves a trade-off between 
the amount of data stored, and the accuracy of the rep- 
resentation of seismicity one can make using that data 
set. From the distribution P(c) shown in Fig. 2, and 
from the average number of incoming links (A:i„) with a 
given choice of c<, we can estimate the error made in 
throwing out weak links. The average correlation c con- 
tribution from all of the « iVnodo = 8858 incoming links 
that are pruned from any earthquake when imposing the 
threshold c< is 

A^nodc / ^ cP[c)dc ~ ^7V„odo4"" , (21) 

where A is a. constant given by the amplitude of P(c), 
and Cmin is the minimum value of c observed in the mea- 
surement of P{c). 

The average correlation contribution from the incom- 
ing links actually represented in the network with that 
choice of c< is 

ihn) / cP{c)dc ~ A{h^)c^mal , (22) 

with Cmax ~ 10^^ for the 2D metric. The relative error 
with a given c< can thus be estimated as 

Errors ---(- ) . (23) 

\i^i7i / '-■max 

With our choice c< — 10* we get an estimate of 
the relative error to be approximately one percent, or 
Error = 0.013 with the fraction the data set stored 
A^iinfcs /A^nodc ~ 0.002. In other words, throwing out 
about 99.8% of the data set we can accurately repre- 
sent the correlation structure of seismicity using a sparse 
network with an estimated error of order one percent. 
Conversely, we are not aware of any quantitative esti- 
mate of the error with particular choices of space time 
windows. 



C. Bench mark test for models of seismicity 

The statistical properties of the network of seismicity 
we find can be used to test various models of seismic- 
ity. A self-organized critical model proposed by Olami 
et al. [S^ exhibits a universal Gutenberg-Richter law for 
earthquakes p9j | , independent of the dissipation parame- 
ter, as well as foreshocks and aftershocks 16]. However, 
no evident self-organized spatial structure corresponding 
to the recurrence of earthquakes on a heterogeneous sys- 
tem of faults exists. For this reason, we believe it is 
unlikely that this model can reproduce the observed net- 
work properties of seismicity. Although no satisfactory 
dynamical model of the self-organization of the Earth's 
crust and resultant seismicity exists at present, a stochas- 
tic branching process, known as the ETAS model |2HI3^ . 
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or its spatially extended version '14"! could be tested by 
constructing a network using our method for particular 
realizations of that process in space, time and magni- 
tude, and comparing with our results. For instance, the 
appearance of the scaling variable ZIO""^™ combining spa- 
tial distances with main shock magnitude could be as- 
certained. Since the distance variable between mother 
daughter pairs in the spatially extended ETAS model is 
chosen from a power law distribution, $(r), independent 
of the parent's magnitude, this model is unlikely to repro- 
duce observed behavior and would have to be modified. 
Conversely, one could also check if our method of con- 
structing networks linking main shocks and aftershocks 



correctly identifies mot her- daughter pairs given by the 
algorithm of the ETAS process or if there might be dif- 
ferences. 
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